A) Movimiento de un Cometa

Existen muchos cometas que orbitan alrededor del Sol con trayectorias elípticas. De acuerdo con las leyes de Kepler, cuando el cometa se encuentra lejos del centro de fuerza, su movimiento es lento, minetras que al acercarse al sol su velocidad es bastante grande. Este es un sistema físico que claramente debería ser solucionado con un método de paso adaptativo: en la región lejana se pueden utilizar pasos grandes mientras que en la región cercana se necesitan pasos muy cortos para tener un error pequeño-.

Utilizando coordenadas cartesianas en el plano de movimiento, la dinámica del cometa bajo la influencia Solar está descrito por las ecuaciones

\begin{align} \frac{d^2x}{dt^2} &= -GM \frac{x}{r^3}; & \frac{d^2y}{dt^2} = -GM \frac{y}{r^3} \end{align}

donde $r = \sqrt{x^2 + y^2}$, G es la constante gravitacional de Newton y M es la masa del Sol.

1.

Defininimos las constantes del sistema en unidades de masas solares, años y unidades astronómicas. En estas unidades las leyes de Kepler nos muestran que la constante de Cavendish toma el valor teórico

\begin{equation} G = 4\pi^2\;. \end{equation}

Luego, reescribimos las dos ecuaciones diferenciales de tal manera que se obtenga un sistema de cuatro ecuaciones diferenciales acopladas de orden uno, como sigue

\begin{align} \frac{dx}{dt} &= v_x; & \frac{dy}{dt} &= v_y;\\ \frac{dv_x}{dt} &= -GM\frac{x}{r^3}; & \frac{dv_y}{dt} &= -GM\frac{y}{r^3} \end{align}

Así, definimos una función que incluya el lado derecho de estas ecuaciones diferenciales

Luego, teniendo en cuenta que el alogoritmo Runge-Kutta de orden 4 es[1]

\begin{align} k_1 &=\Delta x f(x_n,y_n)\,\,,\\ k_2 &=\Delta x f(x_n + \frac{\Delta x}{2},y_n + \frac{1}{2}k_1)\,\,,\nonumber\\ k_3 &=\Delta x f(x_n + \frac{\Delta x}{2},y_n + \frac{1}{2}k_2)\,\,,\nonumber\\ k_4 &=\Delta x f(x_n+\Delta x,y_n + k_3)\,\,,\nonumber\\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2 k_2 + 2 k_3 + k_4) + \mathcal{O}(\Delta x^5)\,\,, \end{align}

definimos una función para desarrollar este algoritmo

Luego, definiendo las condiciones iniciales $t = 0$ y

\begin{array} \\x(0) &= 4 \times 10^9\;km\\ y(0) &= 0\;km\\ vx(0) &= 0 \;m\;s^{-1}\\ vy(0) &= 500 \;m\;s^{-1}\,, \end{array}

que reescalando para que tengan las mismas unidades de las constantes mencionadas previamente, tenemos que estas son

\begin{array} \\x(0) &= 26.738\;ua\\ y(0) &= 0\;ua\\ vx(0) &= 0 \;ua\;yr^{-1}\\ vy(0) &= 0.105 \;ua\;yr^{-1}\,, \end{array}

obtenemos que,

De estos resultados es importante resaltar que el tamaño del paso tomado es importante puesto que dependiendo de este, la simulación puede hacer que el cometa no converja a una órbita cerrada, sino que se escape del sistema, lo cual no representaría correctamente la situación física real, como se puede ver en la siguiente celda

Ahora, para verificar la conservación de cantidades conservadas, definimos dos funciones

Estas cantidades son esencialmente las densidades de momentum angular y energía puesto que, al ser la masa del cometa despreciable respecto a la masa del sol, no tiene sentido físico verificar el comportamiento neto de estas cantidades. Así, tenemos que

2.

Ahora, tenemos el algoritmo de paso adaptativo[1]

\begin{equation} \begin{aligned} k_1 &= \Delta x f(x_n, y_n)\,\,,\\ k_2 &= \Delta x f(x_n + \frac{1}{2} \Delta x, y_n + \frac{1}{2} k_1)\,\,,\\ k_3 &= \Delta x f(x_n + \frac{3}{4} \Delta x, y_n + \frac{3}{4} k_2)\,\,,\\ y_{n+1} &= y_n + \frac{2}{9} k_1 + \frac{1}{3} k_2 + \frac{4}{9} k_3 + \mathcal{O}(\Delta x^4)\,\,\\ k_4 &= \Delta x f(x_n + \Delta x, y_{n+1})\,\,\\ y^*_{n+1} &= y_n + \frac{7}{24} k_1 + \frac{1}{4}k_2 + \frac{1}{3}k_3 + \frac{1}{8} k_4 + \mathcal{O}(\Delta x^3)\,\,. \end{aligned} \end{equation}

Y definiendo

\begin{equation} \Theta = \frac{|\delta y_{n+1}|}{\epsilon}\,\,, \end{equation}

siendo $\delta y_{n+1} = y_{n+1} - y^*_{n+1}\,$ y $\epsilon$ el valor de tolerancia, se determinan las condiciones de actualización del paso como

  1. Considere un paso $\Delta x$ y calcule el valor de $\Theta$.

  2. Si $\Theta > 1$ (error muy grande), entonces

    • defina $\Delta x_\text{new} = \Delta x \left| \frac{1}{\Theta} \right|^{\frac{1}{p}} S$, donde $S$ es un factor de suavizado ($\sim$ $0.9$ ).

    • Descarte la iteración anterior y repitala con $\Delta x_\mathrm{new}$.

  3. Si $\Theta < 1$ (error muy pequeño), entonces

    • defina $\Delta x_\text{new} = \Delta x \left| \frac{1}{\Theta} \right|^{\frac{1}{p}} S$.

    • Acepte la iteración anterior y realice la siguiente iteración con $\Delta x_\text{new}$.

Definimos inicialmente una función de integración Runge-Kutta 3/4 para un paso

Y luego una función para el paso adaptativo

Donde se empleó un ciclo while en vez de un ciclo for puesto que al no tener un paso definido, no se puede establecer un array fijo de valores, sino que estos se obtienen del paso adaptativo. Así, aplicando las condiciones iniciales, una valor de iteración máxima de $2\times 10^4$ y un valor de tolerancia $1\times 1^{-8}$, obtenemos los valores de posición y velocidad actualizados, junto con el array de los respectivos valores de tiempo con uno de los pasos utilizados para hallar cada valor de posición

Donde la relación con la posición de la órbita junto con la actualización del paso se muestra en la siguiente gráfica.

Mostrando que el paso se estaba actualizando adecuadamente.

Por último, procedemos a observar el comportamiento de las cantidades conservadas

Errores porcentuales los cuales son menores que en el método de Runge-Kutta de orden 4, i.e., en el Runge-Kutta 3/4 con paso adaptativo se conservan mejor el momentum y la energía total.

Es importante revisar también si computacionalmente hablando el paso adaptativo es más eficiente. Entonces encontramos que

Vemos que runge kutta adaptativo se demora 10 veces menos que Runge Kutta con paso fijo.

B Orbita Terrestre

Las ecuaciones para el movimeinto de la Tiera alrededor del Sol son

\begin{equation} \begin{split} &\frac{d^2\vec{r}}{dt^2}=-GM\frac{\vec{r}}{r^3} \end{split} \end{equation}

donde G es la constante gravitacional de Newton y M es la masa del Sol. Como es bien conocido, la órbita terrestre no es perfectamente circular. En el punto de máximo acercamiento al Sol (perihelio), la distancia entre los dos cuerpos es de $1,4710 \times 10^{11}$ m y su velocidad lineal es de $3,0287 \times 10^4$ m/s.

1.

Escriba un programa que calcule la orbita terrestre utilizando uno de los métodos simplécticos con un paso temporal de 1hora. Grafique la trayectoria completando varias vueltas alrededor del Sol.

Para hacer la integración se va a utilizar el método de integración PEFRL el cual es simplectico y además de orden 4 en la posición y de orden 3 en la velocidad. Hay una variante de este mismo método que es de orden 4 en la velocidad y orden 3 en la posición. El método está dado por

\begin{equation} \label{PEFRL} \begin{split} &\vec{r}_1=\vec{r}(t)+\vec{v}(t)\zeta dt,\\ &\vec{v}_1=\vec{v}(t)+\frac{1}{2m}\vec{F}(\vec{r}_1)(1-2\lambda)dt,\\ &\vec{r}_2=\vec{r}_1+\vec{v}_1\chi dt,\\ &\vec{v}_2=\vec{v}_1+\frac{1}{m}\vec{F}(\vec{r}_2)\lambda dt,\\ &\vec{r}_3=\vec{r}_2+\vec{v}_2(1-2(\chi+\zeta))dt,\\ &\vec{v}_3=\vec{v}_2+\frac{1}{m}\vec{F}(\vec{r}_3)\lambda dt,\\ &\vec{r}_4=\vec{r}_3+\vec{v}_3\chi dt,\\ &\vec{v}(t+dt)=\vec{v}_3+\frac{1}{2m}\vec{F}(\vec{r}_4)(1-2\lambda)dt,\\ &\vec{r}(t+dt)=\vec{r}_4+\vec{v}(t+dt)\zeta dt,\\ \\ &\zeta=0.1644986515575760E+00,\\ &\lambda=-0.2094333910398989E-01,\\ &\chi=0.1235692651138917E+01. \end{split} \end{equation}

$\vec{v}(t+dt)$ y $\vec{r}(t+dt)$ son la posición y la velocidad en el siguiente paso de tiempo. El artículo con la deducción del método y la demostración de que es simplectico es

Omelyan, I. Mryglod, and R. Folk, Optimized forest–ruth and suzuki-like algorithms for integration of motion in many-body systems, Computer Physics Communications 146, 188 (2002).

Voy a hablar un poco del diseño del programa. El struct Molecule representa cada cuerpo de la simulación. De esta manera la evolución temporal se puede hacer de manera más organizada. Para la posición, velocidad y fuerza se utilizan StaticArrays. Estos Arrays tienen tamaño fijo, pero son los más rápidos que hay. Además las operaciones entre ellos están vectorizadas. Como los vectores velocidad, posición y fuerza siempre tienen 3 componentes estos arrays son ideales para estas cantidades.

Luego se definen los métodos necesarios para actualizar cada cuerpo. La definición de la mayoría de ellos está hecha para actualizar un solo elemento, sin embargo, cuando una función de estas recibe un array y se le coloca un punto, por ejemplo move_v.(body,dt,const2), la operación de actualización se realiza elementwise de forma vectorizada, obteneiendo una gran reducción en overhead de la operación. Por último, se separa la función que hace la evolución temporal de la que hace el paso de tiempo para que cambiar el algoritmo de integración se reduzca simplemente a crear otra función paso y reemplazar el nombre dentro de la propagación. Aprovechando este mismo comportamiento vectorial, se pueden inicializar todas las partículas a la vez como se puede ver a continuación.

2.

Incluya en el programa una función que calcule la energía potencial, la energía cinética y la energía total del sistema en cada paso y grafique sus resultados en los mismos ejes. Los resultados deben mostrar que las energías cinética y potencial cambian visiblemente a lo largo de la trayectoria mientras que la energía total debe permanecer aproximadamente constante.

El programa guarda toda la información en un archivo y solo tiene en memoria la información de cada partícula. Es decir, tiene complejidad espacial O(N), entonces hay que leer la información del archivo y pasarsela a la función que va a calcular la energía en cada paso de tiempo guardado.

3.

Ahora grafique únicamente la energía total para comprobar que existe una pequeña variación en cada una de las orbitas. Sin embargo, el método simpléctico debe retornar al valor inicial después de cada vuelta.

La variación en energía es tan pequeña que el programa tiene problemas para graficar la energía.

4.

Utilice el programa que escribió para verificar si puede describir el movimiento del planeta enano Plutón. Este objeto tiene una orbita mucho mas excentrica que la orbita terrestre, con una distancia al Sol en el perihelio de $4.4368 \times 10^{12}$ m y una velocidad lineal en este punto de $6.1218 \times 10^3$ m/s.

C. Problema Gravitacional de los N-Cuerpos

Las ecuaciones de movimiento de N-partículas moviendose bajo su interacción gravitacional mutua se pueden escribir en la forma \begin{equation} m_i\ddot{x}_i=-Gm_i\sum^N_{j=1,\;i\neq j}\frac{m_j}{|x_{ij}|^2}\hat{x}_{ij}, \end{equation}

donde $x_{ij} = x_i - x_j$ es el vector que apunta desde la partícula j a la partícula i.

1.

Implemente un programa que resuelva el problema diferencial de las ecuaciónes de movimiento utilizando un método simplectico de orden 4. El algoritmo debe ser lo suficientemente general para poder incluír un número arbitrario de partículas y las condiciones iniciales deben ser leídas de un archivo.

El algoritmo del punto anteriror es lo suficientemente general como para hacer la propagación de N moléculas. Sin embargo, todavía falta hacer una función que lea de un archivo las condiciones iniciales e inicialize el sistema. Entonces:

2.

Con el fin de probar el código implementado, utilice los datos iniciales para el sistema SolTierra dados en el archivo sun_earth.dat y grafique el comportamiento del sistema a lo largo de algunos años. Comprube que la orbita no se comporta a la forma de una espiral.

Las funciones diseñadas en el punto anterior ya son generales. Sin embargo hay un detalle. El archivo puede tener datos en unidades distintas a las deseadas en el programa. Es por esto que se hace una función para hacer la conversión de unidades de los datos leidos.

3.

Con el fin de asegurar la convergencia del algoritmo, implmente una rutina dentro del programa que calcule la energía total del sistema de N-partículas. Evalue el comportamiento de la energía a lo largo de la evolución y grafique para comprobar si existe algun cambio significativo. El comportamiento mejora o empeora al modificar el paso de la integración?

Es la misma funcion del ejercicio anterior. Esta ya está generalizada para N cuerpos. Antes con un $dt=0.0005$ la diferencia de energías era de $\Delta E=1.5445040932105986\times 10^{-13}$ %. Un dt más pequeño en general debería mejorar esto, sin embargo, un dt demasiado pequeño abre la puerta para inestabilidades numéricas.

4.

Ahora utilice su código para estudiar la evolución del sistema de 13 estrellas S0 moviendose alrededor del agujero negro supermasivo SgrA*, proporcionado en el archivo S0stars.dat. Verifique también el comportamiento de la energía del sistema de 13+1 partículas. Grafique las orbitas de estas estrellas a lo largo de un periodo de al menos 100 años.

readdlm tiene un problema con los delimitadores, si el archivo tiene espacios iniciales crashea. https://github.com/JuliaLang/julia/issues/44846. Aún así hay un workarround y se pueden cargar los datos sin problemas. La función start_from_file ya tiene implementado el workarround.

También se puede hacer una animación del movimiento de las estrellas. Puede descargar el Gif del repositorio.

5.

Modifique el código para incorporar el método de Verlet para solucionar el problema y resuelva de nuevo para el sistema de las 13 estrellas S0 moviendose alrededor de SgrA*. Mejora el tiempo de computo utilizando este método? Se conserva la energía total del sistema?

Desde el principio se utilizó el método PEFRL que es un método más elaborado que Verlet, sin embargo, involucra muchas más operaciones, Utilizando Verlet y comparando con PEFRL vemos que cambia.

La energía del sistema de 13 estrellas con PEFRL es

Utilizando los mismos parámetros pero con position Verlet tenemos que

Veamos ahora como se comporta velocity Verlet

Vemos que así se gane en un poco en tiempo de ejecución (que no es mucho), para obtener buenos resultados es necesario disminuir el dt de los métodos Verlet, por lo tanto, la mejora en tiempo de ejecución se pierde debido a que toca hacer más pasos de tiempo para obtener el mismo resultado.

Integration method Execution time Energy diference
Position Verlet 0.173431 s 0.06 %
Velocity Verlet 0.228064 s 0.005 %
PEFRL 0.298807 s 3e-8 %

Los tiempos de ejecución puenden ser algo engañosos a primera vista porque incluyen el tiempo de compilación. Para ver únicamente el tiempo de ejecución hacemos un @benchmark

Vemos que realmente no se gana mucho en el tiempo de ejecución pero Verlet pierde mucho en presición.

REFERENCES

[1] Larrañaga, E. Jupyter-notebooks. Astrofísica Computacional.